%[T1Est, bEst, aEst, Res,sd] = MOLLI(data, nlsS)
function [T1Est, bEst, aEst, Res,sd] = MOLLI(data, nlsS)
    
if ( length(data) ~= nlsS.N )
  error('nlsS.N and data must be of equal length!')
end

% Make sure the data come in increasing TI-order
[tVec,order] = sort(nlsS.tVec); 
data = squeeze(data); 
data = data(order);
polluted_pos = zeros(1,length(data));

for ix = 2:length(data)
    if(data(ix)<= data(ix-1))
       polluted_pos(1,ix)  = 1;
       polluted_pos(1,ix-1)  = -1;
    end
end
polluted_pos( tVec<2000 ) = 0;
polluted_pos(end) = 0;
ti = tVec; %( polluted_pos ==0 );
data= data(:);
dat = data; %( polluted_pos ==0 ); 
% x0 = [ max(dat) -2*max(dat) 1400]; 
% [aEst, bEst,T1Est,res]=invdatafit(ti,dat);
% 
% 
% if( sum(abs( polluted_pos ))> 0)
%     [Binary,Pos] = find( polluted_pos==-1 );
%     
%     for ix = 1:length( Pos)
%         ixPos = Pos(ix);
%         Res = zeros(1,2);
%         for i = 1:2
%             ixPos = ixPos + ( i -1);
%         
%             temp_Polluted_pos = polluted_pos;
%             temp_Polluted_pos(1,ixPos) = 0;
% 
%             ti = tVec( temp_Polluted_pos ==0 );
%             dat = data( temp_Polluted_pos ==0 );
%             [aEst, bEst,T1Est,res]=invdatafit(ti,dat);
%             Res(1,i) = res;
%         end
%         %make decision to remove which value
%         if(Res(1,1)>= Res(1,2)  )
%             polluted_pos(1,ixPos) = 0;
%         else
%             polluted_pos(1,ixPos-1) = 0;
%         end
%         %fit fit with the second data
%     end 
% end
% 
% ti = tVec( polluted_pos ==0 );
% dat = data( polluted_pos ==0 );
[aEst, bEst,T1Est,Res]=invdatafit(ti,dat);

%% Find the min of the data

% cal sd
%T1Est = T1Est*(-bEst/aEst - 1);
modelValue = aEst + bEst*exp( - ti/T1Est );
Res = 1/sqrt(length(ti))*norm(1 - modelValue./dat(:));
sd = median( abs( dat - modelValue ))/0.6745 ;
function [aEst, bEst,T1Est,res]=invdatafit(ti,dat)
    
        
    
        
        [Ind, tag] = find(ti<3000);
        aEstHat  = zeros(length(Ind),1);
        bEstHat = zeros(length(Ind),1);
        T1EstHat = zeros(length(Ind),1);
        resHat = zeros(length(Ind),1);
        for ind = 1:length( Ind )
            
            minInd =Ind(ind);
            dataTmp = dat.*[-ones(minInd,1); ones(length(dat) - minInd,1)];
            x0 = [ max(dat) -2*max(dat) 1400];   
            [aEst, bEst,T1Est]= subfitting(ti,dataTmp,x0);
            aEstHat(ind) = aEst;
            bEstHat(ind) = bEst;
            T1EstHat(ind) = T1Est ;%*(-bEstHat(ind)/aEstHat(ind) - 1 );
            modelValue = aEstHat(ind) + bEstHat(ind)*exp( - ti/T1EstHat(ind) );
            resHat(ind) = 1/sqrt(length(dat)) * norm(1 - modelValue./dataTmp);

        end
        [v,ind] = min(resHat);

        aEst = aEstHat(ind(1),1);
        bEst = bEstHat(ind(1),1);
        T1Est = T1EstHat(ind(1),1);
        res = resHat(ind(1),1);
end
function [aEst, bEst,T1Est]=subfitting(ti, data,x0 ) 
        if( length(x0) == 3 )
            try
                x = fminsearch( ...
                @(x)sum( abs( data- (   x(1)*(1+x(2)*exp(-ti/x(3))))).^2 ), ...
            x0,optimset('display','off')); 
            
            catch
                disp('there is an error');
                x=x0 ;
            end
            aEst = x(1);
            bEst = x(2)*x(1);
            T1Est = x(3);
        else
            try
                x = fminsearch( ...
                  @(x)sum( abs( data- (   x(1)*(1-2*exp(-ti/x(2))))).^2 ), ...
                  x0,optimset('display','off')); 
            catch
                  disp('there is an error');
                  x0
                  x = x0; 
                  x
            end 
            aEst = x(1);
            bEst = -2*x(1);
            T1Est = x(2);
        end
end
end
